/*
    Copyright 2007-2012 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Implementation of the ir_stage class. \ingroup mcrx

#include "config.h"
#include "mcrx-stage.h"
#include "mcrx-stage-impl.h"
#include "xfer_impl.h"
#include "create_grid.h"
#include "arepo_grid.h"
#include "arepo_emission.h"

template <template <typename> class grid_type> 
void mcrx::ir_stage<grid_type>::calculate_dust_SED()
{
  CCfits::ExtHDU& intensity_hdu = open_HDU(*this->output_file_, "INTENSITY");
  CCfits::ExtHDU& deposition_hdu = open_HDU(*this->output_file_, "DEPOSITION");

  // get the grain model
  std::vector<boost::shared_ptr<grain_model<polychromatic_scatterer_policy, 
    mcrx_rng_policy> > > models;
  models.push_back
    (load_grain_model<polychromatic_scatterer_policy, mcrx_rng_policy>(this->p_, this->m_.units));
  models[0]->resample(ilambda_, elambda_);
  // trigger generation of interpolation table
  //models[0]->calculate_SED(lambda_,0);

  mcrx_terminator term(this->m_);
  const int n_threads = this->p_.defined("n_threads")?
    this->p_.getValue("n_threads", int ()): 1;

  this->emi_->calculate_SED(intensity_hdu,
			    models, term, n_threads,
			    0, 
			    false, true);
}

template <template <typename> class grid_type> 
void mcrx::ir_stage<grid_type>::basic_setup ()
{
  // read original and intensity wavelength vectors
  CCfits::ExtHDU& lambda_hdu = open_HDU (*this->output_file_, "LAMBDA");
  read(lambda_hdu.column("lambda"), elambda_);
  read(lambda_hdu.column("lambda_intensity"), ilambda_);
  int n_lambda_int;
  lambda_hdu.readKey("NLAMBDA_INTENSITY", n_lambda_int);
  ilambda_.resizeAndPreserve(n_lambda_int);

  // Create ir_grid object which calculates the dust emission in the
  // grid cells
  cout << "Setting up emission" << endl;

  // create density generator
  uniform_density_generator dg(this->p_);

  // create ir_grid as the emission object
  const int decomp_level = this->p_.getValue("domain_level", int(1), true);
  T_grid* dummy=0;
  this->emi_ = create_grid(*this->input_file_, dg, elambda_, 
			   decomp_level, dummy);

  // Build cameras based on -PARAMETERS HDUs in output file
  cout << "Setting up emergence" << endl;
  this->eme_.reset(new T_emergence (*this->output_file_));

  // if we are doing kinematic calculation, set the wavelength scale
  // of the cameras
  if(this->p_.getValue("use_kinematics", false, true)) {
    typename T_emergence::iterator stop= this->eme_->end();
    for (typename T_emergence::iterator c = this->eme_->begin();
	 c != stop; ++c)
      (*c)->set_dopplerscale(elambda_);
  }

  // tell cameras to allocate
  typename T_emergence::iterator stop= this->eme_->end();
  for (typename T_emergence::iterator c = this->eme_->begin();
       c != stop; ++c) {
    (*c)->allocate(elambda_);
    assert((*c)->get_image().size()>0);
  }

}

template <template <typename> class grid_type> 
void mcrx::ir_stage<grid_type>::setup_objects ()
{
  basic_setup();
  calculate_dust_SED();
}



/** Load IR stage data from a previously completed file. This is not
    so efficient because we don't save the emission SEDs of the grid
    cells, so the temperature distribution has to be recalculated
    either way. */
template <template <typename> class grid_type> 
void mcrx::ir_stage<grid_type>::load_file ()
{
  basic_setup();
  calculate_dust_SED();

  if (!this->m_.terminator()) {
    // can only load images, sed is not saved
    this->eme_->load_images(*this->output_file_, "-"+stage_name());
  }
}


template <template <typename> class grid_type> 
void mcrx::ir_stage<grid_type>::load_dump (binifstream& dump_file)
{
  basic_setup();

  if (!this->m_.terminator()) {
    // Load dust temp distribution, cells remaining and the current images
    this->emi_->load_dump(dump_file);
    this->eme_->load_dump(dump_file);
    // finish temp calculation
    calculate_dust_SED();
  }
}


template <template <typename> class grid_type> 
void mcrx::ir_stage<grid_type>::save_file ()
{
  // Get normalization and save images
  //const T_float normalization = this->emi_->total_emission_weight()/this->n_rays_;
  const T_float normalization = 1.0;
  this->eme_->write_images(*this->output_file_, normalization, this->m_.units,
		     "-"+stage_ID(),
		     this->p_.defined("compress_images")&&
		     this->p_.getValue("compress_images", bool ()),
		     this->p_.defined("write_images_parallel")&&
		     this->p_.getValue("write_images_parallel", bool ()));

}


/** Dump IR stage info. We might have been interrupted during the temp
    calculation, or we might have been interrupted during shooting. In
    either case, we should save the cell SEDs that have been
    calculated so we don't have to do it again. */
template <template <typename> class grid_type> 
void mcrx::ir_stage<grid_type>::save_dump (binofstream& dump_file) const
{
  this->emi_->write_dump(dump_file);
  this->eme_->write_dump(dump_file);
}


template <template <typename> class grid_type> 
bool mcrx::ir_stage<grid_type>::shoot ()
{
  return this->shoot_nonscatter ();
}

template <template <typename> class grid_type> 
void mcrx::ir_stage<grid_type>::operator() ()
{
  // this is here so the function is instantiated
  this->run_stage();
}



// prevent instantiation of the xfer template in each stage. instead
// we do this in one place only. (The dummy grid version of xfer in
// implicitly instantiated because we can't explicitly instantiate it
// since not all methods compile with the dummy grid.)
extern template class 
mcrx::xfer<mcrx::dust_model<mcrx::scatterer<mcrx::polychromatic_scatterer_policy,
					    mcrx::mcrx_rng_policy>,
			    mcrx::cumulative_sampling,
			    mcrx::mcrx_rng_policy>,
	   mcrx::T_full_sed_adaptive_grid>;
#ifdef WITH_AREPO
extern template class 
mcrx::xfer<mcrx::dust_model<mcrx::scatterer<mcrx::polychromatic_scatterer_policy,
					    mcrx::mcrx_rng_policy>,
			    mcrx::cumulative_sampling,
			    mcrx::mcrx_rng_policy>,
	   mcrx::T_full_sed_arepo_grid>;
#endif

// explicit instantiations
template class mcrx::ir_stage<mcrx::adaptive_grid>;

#ifdef WITH_AREPO
template class mcrx::ir_stage<mcrx::arepo_grid>;
#endif
